# "Media's Influence on LGBTQ Support Across Africa" by Stephen Winkler

# Code to create marginal effects plot for the interaction of KOF with media consumption:
#   - Figure 4 in main text

# only need to plot the logit model.
# substantive results are stable across OLS & Multilevel models

rm(list=ls()) # clear environment
source("Rprofile.R") # setup Rprofile (load packages, etc.)

# Load data
myData <- readRDS("data/afrobarometer.rds")

# Convert Country into Dummy & Rename
dummies <- predict(dummyVars(~ country, data = myData, levelsOnly = TRUE), newdata = myData)
colnames(dummies) <- gsub("country", "", colnames(dummies)) #remove country from column names
colnames(dummies) <- gsub(" ", "", colnames(dummies)) #remove spaces
colnames(dummies) <- gsub("'", "", colnames(dummies)) # remove unusual characters

# Join country dummies back to main dataset
myData <- cbind(myData, dummies)
myData <- dplyr::select(myData, -Algeria, -Sudan, -Egypt, -Zambia, -Zimbabwe) #remove NA countries
                            #NOTE: I have data for Zambia & Zimbabwe but they return NA coefficients in the models
                            # this means they are likely missig data for the interaction (?)

myData <- filter(myData, country != "Algeria")
myData <- filter(myData, country != "Sudan")
myData <- filter(myData, country != "Egypt")

# convert all variables to appropriate class 
# can't do it in the model because using language like 'as.numeric' 
# will throw off the cplot() function when I make marginsplots
myData$radio <- as.numeric(myData$radio)
myData$tv <- as.numeric(myData$tv)
myData$newspaper <- as.numeric(myData$newspaper)
myData$internet <- as.numeric(myData$internet)
myData$social_media <- as.numeric(myData$social_media)
myData$media_aggregate <- as.numeric(myData$media_aggregate)
myData$media_noradio <- as.numeric(myData$media_noradio)
myData$media_notv <- as.numeric(myData$media_notv)
myData$media_nopaper <- as.numeric(myData$media_nopaper)
myData$media_nointernet <- as.numeric(myData$media_nointernet)
myData$media_nosocialmedia <- as.numeric(myData$media_nosocialmedia)
myData$education <- as.numeric(myData$education) #if not, then get error with cplot
myData$religiosity <- as.numeric(myData$religiosity)
myData$age <- as.character(myData$age)
myData$age <- as.numeric(myData$age)
myData$income_water <- as.numeric(myData$income_water)
myData$country <- as.factor(myData$country)
myData$female <- as.numeric(myData$female)
myData$urban <- as.numeric(myData$urban)

############################################
# Robustness Check: KOF Interaction
############################################
# Rerun main models but with an interaction for globalization

# Radio Logit
model_radio <-  # model 
  sexuality2 ~
  radio + radio:KOFSoGI +
  media_noradio+ #consumption of other mediums
  tolerance_noLGBT +
  female + education + religiosity + age + income_water + urban +
  factor(country)

mdata_radio <- extractdata(model_radio, myData, na.rm=TRUE) #only extract data used in model
logit.radio <- glm(model_radio, data = mdata_radio, family = binomial)
# Code below uses the margins package to extract the effect of medium across 
# different levels of KOF. This code takes a long time to run so do not rerun it unless necessary.
# results are already saved as an RDS and are read in later in the code. 
#cplot.radio.logit <- cplot(logit.radio, x ="KOFSoGI", dx = "radio", what = "effect", se.type = "shade", draw = FALSE)
#saveRDS(cplot.radio.logit, file = project_path("analysis/results", "cplot.radio.logit.rds"))

# TV Logit
model_tv <- 
  sexuality2 ~
  tv + tv:KOFSoGI +
  media_notv+ #consumption of other mediums
  tolerance_noLGBT +
  female + education + religiosity + age + income_water + urban +
  factor(country)

mdata_tv <- extractdata(model_tv, myData, na.rm=TRUE) #only extract data used in model
logit.tv <- glm(model_tv, data = mdata_tv, family = binomial)
# Code below uses the margins package to extract the effect of medium across 
# different levels of KOF. This code takes a long time to run so do not rerun it unless necessary.
# results are already saved as an RDS and are read in later in the code.
#cplot.tv.logit <- cplot(logit.tv, x ="KOFSoGI", dx = "tv", what = "effect", se.type = "shade", draw = FALSE)
#saveRDS(cplot.tv.logit, file = project_path("analysis/results", "cplot.tv.logit.rds"))

# Newspaper Logit
model_newspaper <- 
  sexuality2 ~
  newspaper + newspaper:KOFSoGI +
  media_nopaper + #consumption of other mediums
  tolerance_noLGBT +
  female + education + religiosity + age + income_water + urban +
  factor(country)

mdata_newspaper <- extractdata(model_newspaper, myData, na.rm=TRUE) #only extract data used in model
logit.newspaper <- glm(model_newspaper, data = mdata_newspaper, family = binomial)
# Code below uses the margins package to extract the effect of medium across 
# different levels of KOF. This code takes a long time to run so do not rerun it unless necessary.
# results are already saved as an RDS and are read in later in the code.
#cplot.newspaper.logit <- cplot(logit.newspaper, x ="KOFSoGI", dx = "newspaper", what = "effect", se.type = "shade", draw = FALSE)
#saveRDS(cplot.newspaper.logit, file = project_path("analysis/results", "cplot.newspaper.logit.rds"))

#Internet Logit
model_internet <- 
  sexuality2 ~
  internet + internet:KOFSoGI +
  media_nointernet+ #consumption of other mediums
  tolerance_noLGBT +
  female + education + religiosity + age + income_water + urban +
  factor(country)

mdata_internet <- extractdata(model_internet, myData, na.rm=TRUE) #only extract data used in model
logit.internet <- glm(model_internet, data = mdata_internet, family = binomial)
# Code below uses the margins package to extract the effect of medium across 
# different levels of KOF. This code takes a long time to run so do not rerun it unless necessary.
# results are already saved as an RDS and are read in later in the code.
#cplot.internet.logit <- cplot(logit.internet, x ="KOFSoGI", dx = "internet", what = "effect", se.type = "shade", draw = FALSE)
#saveRDS(cplot.internet.logit, file = project_path("analysis/results", "cplot.internet.logit.rds"))

# Social Media Logit
model_social_media <- 
  sexuality2 ~
  social_media + social_media:KOFSoGI +
  media_nosocialmedia + #consumption of other mediums
  tolerance_noLGBT +
  female + education + religiosity + age + income_water + urban +
  factor(country)

mdata_social_media <- extractdata(model_social_media, myData, na.rm=TRUE) #only extract data used in model
logit.social_media <- glm(model_social_media, data = mdata_social_media, family = binomial)
# Code below uses the margins package to extract the effect of medium across 
# different levels of KOF. This code takes a long time to run so do not rerun it unless necessary.
# results are already saved as an RDS and are read in later in the code.
#cplot.social_media.logit <- cplot(logit.social_media, x ="KOFSoGI", dx = "social_media", what = "effect", se.type = "shade", draw = FALSE)
#saveRDS(cplot.social_media.logit, file = project_path("analysis/results", "cplot.social_media.logit.rds"))

# combine all the logit cplots together
cplot.radio.logit <- readRDS("analysis/results/cplot.radio.logit.rds")
cplot.tv.logit <- readRDS("analysis/results/cplot.tv.logit.rds")
cplot.newspaper.logit <- readRDS("analysis/results/cplot.newspaper.logit.rds")
cplot.internet.logit <- readRDS("analysis/results/cplot.internet.logit.rds")
cplot.social_media.logit <- readRDS("analysis/results/cplot.social_media.logit.rds")

cplot.all.logit <- bind_rows(cplot.radio.logit, cplot.tv.logit, cplot.newspaper.logit, cplot.internet.logit, cplot.social_media.logit)
cplot.all.logit %<>% # change names to match other figures in the paper
  mutate(factor = recode_factor(factor, "radio" = "Radio",
                                "tv" = "TV",
                                "newspaper" = "Newspaper",
                                "internet" = "Internet",
                                "social_media" = "Social media"))

cplot.all.logit$factor <- factor(cplot.all.logit$factor,
                                 levels = c("Radio", "TV", "Newspaper",
                                            "Internet", "Social media"
                                    )) 

# Plot 
# Figure 4 in Main Text
ame_kof_logit_plot <-
  ggplot(cplot.all.logit, aes(x = xvals)) +
  geom_ribbon(aes(ymin = lower, ymax = upper, group = factor), fill = "gray80") +
  geom_line(aes(y = yvals, 
            group = factor
            #color = factor
            )) +
  scale_x_continuous(breaks = c(40, 70), labels = c("Low", "High")) +
  #geom_hline(yintercept=c(0), linetype="dotted") + 
  xlab("level of social globalization (KOF)") +
  ylab("average marginal effect of medium") +
  labs(caption = "'Low' represents lowest KOF score in the data (35) 
                  'High' represents highest KOF score in the data (73)") + 
  #ggtitle("Average marginal effect of media on LGBTQ support across levels of social globalization") +
  theme_bw() + theme(#panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     #panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5),
                     axis.ticks = element_blank()) +
                     #axis.text.x = element_text(angle=45)
                     #legend.position = c(.9,.5)
  #facet_wrap(~ factor)
  facet_grid(.~factor)
ame_kof_logit_plot
ggsave(filename = "figures/ame_kof_logit.pdf", #Save it to directory
       plot = ame_kof_logit_plot,
       width = 8, height = 6)
